← LAB ↑ TOC ↑

Name, Vorname:
Hofmann Fabian
Richter Daniel

LAB 2: Filterung und DFT als Näherung für spektrale Analyse

Inhalt

0. Allgemeines

1. FIR und IIR Filter

2. Schätzung der Fourierreihe durch FFT

3. Schätzung des Fourierspektrums von Sinustönen mit überlagertem Rauschen durch FFT

4. Fourierspektrum eines Rechteckpulses

0. Allgemeine Hinweise

Nach dem Praktikumsversuch exportieren Sie das Notebook mit Textantworten, Codezellen und Plots als HTML (File -> Export Notebook As ... -> Export Notebook to HTML) und reichen es in Moodle ein.</div>


1. FIR und IIR Filter

In diesem Abschnitt analysieren Sie die beiden Filter aus der folgenden Abbildung:

  
Filter 1  Filter 2
VORBEREITUNG:

1.1 Simulation

Für die Simulation müssen wir zuerst die Filter als Koeffizientenvektoren a und b darstellen (siehe folgendes Bild und Abschnitt 3 von LAB 1).

Als "Codesteinbruch" können Sie wieder das Notebook 02_LTF/LTF-Filter_properties.ipynb verwenden.

Alternativ verwenden Sie pyfda: Im Tab "b,a" importieren Sie Koeffizienten aus einem CSV-File oder (nach Auswahl von "Clipboard" in den Einstellungen) direkt aus der Zwischenablage. Die Koeffizienten b des nicht-rekursiven Teils geben Sie hierfür getrennt durch Kommata an, optional in einer zweiten Zeile die Koeffizienten a des rekursiven Teils.

SIMULATION:

1.1.1 P/N Diagramm

1.1.2 Impulsantwort

Der Befehl h,t = dsp.impz(b,a) berechnet die Impulsantwort, die Sie dann plotten können, am Besten als stem-Plot.

Warum ist beim IIR-Filter jeder zweite Impuls Null?

1.1.3 Betrags- und Phasengang sowie Gruppenlaufzeit

Aus dem komplexen Frequenzgang omega, H = sig.freqz(b,a) ermitteln Sie mit np.abs() und np.angle() Betrags- und Phasengang. Defaultmäßig werden 512 Frequenzpunkte zwischen 0 und $f_S/2$ bestimmt und als normierte Kreisfrequenz zurückgegeben.

Die Gruppenlaufzeit ermitteln Sie mit w, H = sig.group_delay((b,a), omega). Mit omega übergeben Sie die Frequenzen, an denen die Gruppenlaufzeit berechnet werden soll (z.B. die, die Sie aus der Berechnung des komplexen Frequenzgangs erhalten haben).

Der Plot der Gruppenlaufzeit sieht u.U. etwas seltsam aus, mit set_ylim([y_min, y_max]) passen Sie die Grenzen an.

1.2 Anhören

Im folgenden filtern wir Sprach- oder Rauschsignale mit unseren beiden Filtern und hören und schauen uns das Resultat an. Eine FIR-Filterung könnten Sie wieder mit convolve(x,h) durchführen (nur für eindimensionale Arrays). Warum funktioniert bei IIR-Filtern die convolve-Methode nicht?

Für IIR und FIR-Filter können Sie die Routine y = sig.lfilter(b,a,x) verwenden (funktioniert auch mit zweidimensionalen Arrays).

SIMULATION:

2. Fourierreihe einer Rechteckpulsfolge

In diesem Versuchsteil bestimmen wir die Fourierkoeffizienten einer rechteckförmigen zeitkontinuierlichen Pulsfolge mit $T_1 = 1\,$ms. Der Duty Cycle $\alpha = 0.25$ und Amplitude $A=1000\,$mV, die Pulsbreite ist also $\Delta T = T_1/4$ (siehe nächster Plot).

2.1 Berechnung der Koeffizienten

Ein mit $T_1$ periodisches Signal $y(t)$ lässt sich mit Hilfe der Fourierreihe für reellwertige Signale (Index $k \in \mathbb{N}_0$) darstellen:

$$ y(t)= {a_0} + 2\sum_{k=1}^\infty a_k \cos 2 \pi k f_1 t + b_k \sin 2 \pi k f_1 t \; k \in \mathbb{N} $$

Zunächst berechnen wir die Koeffizienten $a_k$, $b_k$ für $k = 0, 1, \ldots, 11$ mit Hilfe der Formel für die Fourierreihe von reellwertigen Zeitsignalen. Da das Signal achsensymmetrisch ist, sind die imaginärwertigen Koeffizienten $b_k=0$.

$$ \begin{align}{a}_k &= \frac 1 {T_1} \int_{-T_1/2}^{T_1/2}y(t)\cos(2 k \pi f_1 t) \, \text{d}t = \frac {A} {T_1} \int_{-\alpha T_1/2}^{\alpha T_1/2}\cos(2 k \pi f_1 t) \, \text{d}t \\ &= \left. \frac {A} {T_1 \cdot 2 k \pi f_1} \sin(2 k \pi f_1 t)\, \right|_{-\alpha T_1/2}^{\alpha T_1/2} = {\frac {A} { k \pi} \sin ( k \alpha \pi) }= { {A \alpha}\, \mathrm{si} ( k \alpha \pi) } \end{align}$$
SIMULATION:

Schreiben Sie ein Skript, das die ersten 11 Koeffizienten der Rechteckfunktion exakt berechnet. Welchen Vorteil hat es, die sinc(x) Funktion zu verwenden anstelle von sin(x)/x? Aber Achtung: In Numpy (und in Matlab) ist die sinc-Funktion definiert als $\mathrm{sinc}(x) = sin(\pi x)/\pi x$, das $\pi$ im Argument müssen Sie daher weglassen. Drucken Sie die Ergebnisse übersichtlich in eine Tabelle aus (s.u.)

2.1.1 Formatierte Ausdrucke in Python

Einen formatierten Ausdruck erhalten Sie z.B. mit print("\n{0:7.2f} | ".format(Y[i]), end="") (der Teil in der geschweiften Klammer wird ersetzt durch Y[i] und formatiert mit insgesamt 7 Stellen, 2 Nachkommastellen, keinen Zeilenumbruch). Eingebettet in eine for Schleife bekommen Sie so schnell eine übersichtliche Tabelle.

2.2 Abschätzung des Spektrums durch FFT

Um das Spektrum des zeitkontinuierlichen Signals mit einer FFT abschätzen zu können, muss es zunächst abgetastet werden und zwar über eine ganzzahlige Anzahl Perioden $T_1$ (da es ja periodisch ist), wir starten zunächst mit einer Periode, $T_{mess1} = T_1$. Für eine effiziente Berechnung sollen $N_1 = 2^4 = 16$ Samples genommen werden.

VORBEREITUNG:
SIMULATION:

Zunächst einmal benötigen Sie einen passenden Zeitvektor z.B. mit t1 = np.arange(N_1)/f_S1 oder durch Abtasten des "zeitkontinuierlichen" Zeitvektors t (s.u.). Das abgetastete Signal $y_1[n]$ erzeugen Sie:

Egal für welche Variante Sie sich entschieden haben, plotten Sie das abgetastete Signal mit dem folgenden Skript.


2.2.1 FFT, erster Ansatz: $f_{S1} = 16 \text{ kHz}$, $N_1 = 16$

VORBEREITUNG:
SIMULATION:

In numpy berechnen wir die FFT mit Y = np.fft.fft(y, N) (die FFT liefert komplexe Werte ...). Ohne den Parameter N wird die FFT über das gesamte Signal berechnet mit der Anzahl der Datenpunkte. Nutzen Sie das Notebook 03_DFT/DFT-Skalierung.ipynb zur korrekten Berechnung und Skalierung des Signals. Mit f = np.fft.fftfreq(N, T_S) erzeugen Sie eine passenden Frequenzvektor.

Nach der SIMULATION:

2.2.2 FFT mit Erhöhung der Abtastfrequenz auf $f_{S2} = 64 \text{ kHz}$ und $N_2$ Samples

Die Abtastfrequenz wird jetzt vervierfacht auf den Wert $f_{S2} = 64 \text{ kHz}$.

VORBEREITUNG:
SIMULATION:

2.2.3 FFT mit verlängertem Messfenster ($T_{Mess3} = 4 \text{ ms}$) und $f_{S3} = 16 \text{ kHz}$

Die Abtastfrequenz ist wieder wie zu Beginn $f_{S3} = 16 \text{ kHz}$, jetzt soll getestet werden wie sich eine Verlängerung des Messfensters um den Faktor 4 auswirkt.

VORBEREITUNG:
SIMULATION:

Da Ihr Messfenster jetzt mehrere Perioden des Signals umfasst, müssen Sie Ihr Signal erzeugen entweder

Plotten Sie Ihr Zeitsignal, damit Sie sicher sind dass Ihr Signal so aussieht wie Sie sich das vorstellen.

2.2.4 FFT über nicht-ganzzahlige Anzahl von Perioden

Das Signal wird wieder mit $f_{S4} = 16\text{ kHz}$ abgetastet, allerdings ist das Messfenster jetzt nur noch $60 \, T_S$ lang.

VORBEREITUNG:

Signal und Zeitvektor erzeugen Sie am einfachsten durch Verkürzen der Signale aus dem letzten Abschnitt, also mit y4 = y3[:60] oder alternativ y4 = y3[:-4] (negative Indizes zählen vom Ende aus). Der Index, der das Ende markiert ("60" bzw. "-4") definiert dabei den ersten Wert, der nicht mit ausgegeben wird. Das ist z.B. bei arange(0,10,1) genauso und ist eine Folge des 0-based indexing (der erste Wert eines Arrays wird mit 0 angesprochen).

3. Schätzung des Fourierspektrums von Sinustönen mit überlagertem Rauschen durch FFT

In diesem Versuchsteil versuchen wir mit Hilfe der FFT Amplituden und Frequenzen von zwei Sinustönen zu schätzen mit $f_1 = 990\text{ Hz}$ und $A_1 = 100\text{ mV}$ sowie $f_2 = 1010\text{ Hz}$ und $A_2 = 2\text{ mV}$, überlagert ist eine gleichverteilte Rauschstörung im Bereich $\pm 1 \text{ mV}$.

Tipp: Diesen Versuchsteil können Sie auch mit pyfda durchführen, im Tab y[n] wählen Sie dazu ein Sinussignal mit überlagertem Rauschen aus und schauen sich das Spektrum im Untertab "Frequency" an. Enablen Sie "Stimulus X" und deaktivieren Sie "Response Y".

3.1 Abtastfrequenz frei wählbar

Die Anzahl der Datenpunkte $N$ für die FFT und die Abtastfrequenz $f_S$ können zunächst frei gewählt werden.

VORBEREITUNG:


SIMULATION:

3.2 Feste Abtastfrequenz $f_S = 4 \text{ kHz}$

Die Abtastfrequenz ist jetzt fest eingestellt auf $f_S = 4 \text{ kHz}$, die Anzahl der Datenpunkte $N$ für die FFT soll eine Zweierpotenz sein, $N = 2^L$, für eine möglichst effiziente Berechnung der FFT.

Auch diesen Unterpunkt können Sie mit pyfda bearbeiten.

VORBEREITUNG:
SIMULATION:

4. Schätzung des Fourierspektrums eines rechteckförmigen Impulses durch FFT

In diesem Versuchsteil schätzen wir das Fourierintegral des rechteckförmigen Impulses $y(t)$ mit $\Delta T = 250 \,\mu$s und $A = 1000\,\text{mV}$ (s. nächster Plot) mit Hilfe der FFT.

4.1 Berechnung der Fouriertransformierten

VORBEREITUNG:

Das zeitkontinuierliche Signal $y(t)$ hat eine endliche zeitliche Ausdehnung und endliche Energie, daher kann sein Amplitudenspektrum mit dem Fourierintegral beschrieben werden:

$$Y\left( f \right) = \int_{-\infty}^{\infty}y(t)\text{e}^{-\text{j} 2 \pi f t} \text {d}t$$

Der Impuls ist symmetrisch zur y-Achse was die Berechnung deutlich vereinfacht:

$$ \begin{align} Y\left( f \right) &= A \int_{-\Delta T/2}^{\Delta T/2}\text{e}^{-\text{j} 2 \pi f t} \text {d}t = \left. \frac A{-\text{j} 2 \pi f} \text{e}^{-\text{j} 2 \pi f t}\right|_{-\Delta T/2} ^{\Delta T/2} = A\frac{\text{e}^{\text{j} \pi f \Delta T}- \text{e}^{-\text{j} \pi f \Delta T}}{2 \pi \text{j}f} = A\frac{\sin \pi f \Delta T}{\pi f} \\ &= A\Delta T \frac{\sin \pi f \Delta T}{\pi f \Delta T} = A\Delta T \text{sinc} f \Delta T \text{ mit } \text{sinc} f = \frac{\sin \pi f}{\pi f} \end{align}$$

Das Ergebnis ist die Amplitudenspektrumsdichte in V/Hz.

Berechnet man die FFT anstatt des Fourierintegrals, überstreicht jeder Frequenzpunkt die Bandbreite $\Delta f = f_S/N$. Möchte man die FFT eines Energiesignals so skalieren dass man (in etwa) identische Zahlenwerte erhält wie bei der spektralen Amplitudendichte des ursprünglichen zeitkontinuierlichen Signals, so muss man mit $1/\Delta f = N/f_S$ multiplizieren. Da die FFT bereits mit $1/N$ skaliert ist, hebt sich der Faktor $N$ auf. Man muss also nur durch $f_S$ dividieren.

Für unendlich ausgedehnte Leistungssignale (periodische Signale, stationäre und nicht-stationäre Prozesse) konvergiert das Integral nicht, man nimmt stattdessen die Fourierreihe (periodische Signale) oder die spektrale Leistungsdichte, die über die Autokorrelationsfunktion berechnet wird (aber nicht hier und jetzt ...).

4.2 Schätzung der Fouriertransformierten mit Hilfe der FFT

Im folgenden soll die Fouriertransformierte des Impulses mit Hilfe einer FFT abgeschätzt werden, dazu wird der Impuls mit $L = 16$ Samples abgetastet.

VORBEREITUNG:

(c) 2016 - 2021 Prof. Dr. Christian Münker

This jupyter notebook is part of a collection of notebooks on various topics of Digital Signal Processing. The latest version can be found at https://github.com/chipmuenk/dsp.

This notebook is provided as Open Educational Resource. Feel free to use it for your own purposes. The text is licensed under Creative Commons Attribution 4.0, the code of the IPython examples under the MIT license. Please attribute the work as follows: Christian Münker, Digital Signal Processing - Vorlesungsunterlagen mit Simulationsbeispielen, 2020.